function hfile = srf_Reduce(hfile)
% SRF::Reduce  - reduce regular (icosahedron) SRF by factor 4
%
% FORMAT:       srf.Reduce;
%
% No Input/Output fields.

% Version:  v0.7f
% Build:    8110521
% Date:     Nov-05 2008, 9:00 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% check arguments
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'srf')
    error( ...
        'BVQXfile:BadArguments', ...
        'Invalid call to %s.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);

% check for regularity
pfac = log2((bc.NrOfVertices - 2) / 10) / 2;
ffac = log2(bc.NrOfTriangles / 20) / 2;
if pfac ~= fix(pfac) || ...
    ffac ~= fix(ffac) || ...
    pfac ~= ffac || ...
    pfac < 1
    error( ...
        'BVQXfile:BadObject', ...
        'Irregular SRF mesh specified. Can''t reduce by factor 4.' ...
    );
end

% build new size
pnum = (2 ^ (pfac * 2 - 2)) * 10 + 2;
fnum = (2 ^ (pfac * 2 - 2)) * 20;

% set new options
bc.NrOfVertices = pnum;
bc.NrOfTriangles = fnum;
bc.VertexCoordinate(pnum+1:end, :) = [];
bc.VertexNormal(pnum+1:end, :) = [];
bc.VertexColor(pnum+1:end, :) = [];
bc.AutoLinkedMTC = '';

% building new connections
cnx = cell(pnum, 2);
ocn = bc.Neighbors(:, 2);
rfc = 0;
for c = 1:pnum
    oc = ocn{c};
    nc = [];
    ii = [];
    for vc = oc(:)'
        iv = intersect(oc, ocn{vc});
        for ic = iv(:)'
            ii = union(ii, intersect(ocn{vc}, ocn{ic}));
        end
        ii = setdiff(ocn{vc}, union(ii, iv));
        nc = [nc ii(:)'];
    end
    if any(nc > pnum)
        error( ...
            'BVQXtools:MathError', ...
            'Invalid point in new connection list entry.' ...
        );
    end
    cnx{c, 1} = length(nc);
    cnx{c, 2} = nc;
    rfc = rfc + length(nc);
end
bc.Neighbors = cnx;

% building new faces
cnx = cnx(:, 2);
nfs = zeros(rfc, 3);
rfc = 1;
for c = 1:pnum
    fcp = [cnx{c} cnx{c}(1)];
    for cc = 1:length(fcp)-1
        if all(fcp(cc:cc+1) > c)
            nfs(rfc, :) = [c fcp(cc:cc+1)];
        elseif fcp(cc) < fcp(cc+1)
            nfs(rfc, :) = [fcp(cc:cc+1) c];
        else
            nfs(rfc, :) = [fcp(cc+1) c fcp(cc)];
        end
        rfc = rfc + 1;
    end
end
bc.TriangleVertex = unique(nfs, 'rows');
if size(bc.TriangleVertex, 1) ~= fnum
    error( ...
        'BVQXtools:MathError', ...
        'Failure reconnecting vertices to faces.' ...
    );
end

% set back
bvqxfile_setcont(hfile.L, bc);

% recalc normals
srf_RecalcNormals(hfile);
